arXi v: 0907. 203 lvl  [math.NA]  13  Jul  2009 


Synthetic  Aperture  Sonar  Imaging  via  One-Way  Wave 

Equations 

Quyen  Huynh*  Kazufumi  Itct 

May  21,  2010 


Abstract 

We  develop  an  efficient  algorithm  for  Synthetic  Aperture  Sonar  imaging  based  on 
the  one-way  wave  equations.  The  algorithm  utilizes  the  operator-splitting  method 
to  integrate  the  one-way  wave  equations.  The  well-posedness  of  the  one-way  wave 
equations  and  the  proposed  algorithm  is  shown.  A  computational  result  against  real 
field  data  is  reported  and  the  resulting  image  is  enhanced  by  the  BV-like  regularization. 


1  Introduction 

In  this  paper  we  discuss  the  migration  method  based  on  one-way  wave  equations  DEI 
for  Synthetic  Aperture  Sonar  (SAS)  imaging  [3],  The  one-way  wave  equation  integrates 
the  data  within  a  given  angle  and  minimizes  the  undesirable  effects  of  unwanted  reflections. 
Efficient  and  stable  integration  methods  of  the  one-way  wave  equation  based  on  the  operator 
splitting  method  are  used  to  develop  a  fully  discretized  algorithm.  The  stability  analysis 
and  the  required  operation  count  of  the  proposed  algorithm  are  given.  We  test  the  proposed 
method  for  real  field  data  and  report  our  SAS  imaging  results.  We  also  discuss  the  image 
enhancement  method  for  the  resulting  images,  based  on  BV-like  regularization  technique  [5J. 

In  side-scan  (side-looking)  sonar  systems  a  platform  containing  a  moderately  large  real 
aperture  antenna  travels  along  a  rectilinear  path  in  the  along  track  direction  and  periodically 
transmits  a  pulse  at  an  angle  that  is  perpendicular  to  the  platform  path.  These  systems 
produce  strip-map  (SAS)  images  .  A  strip-map  image  is  built  up  as  follows;  the  imaging 
system  operates  such  that  the  echoes  from  the  current  pulse  are  received  before  the  next 
pulse  is  transmitted.  As  these  echoes  are  received  they  are  demodulated,  pulse  compressed, 
and  detected  (only  the  magnitude  information  is  retained).  Each  detected  pulse  produces  a 
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range  line  of  the  real  aperture  image.  As  the  platform  moves  these  range  lines  are  displayed 
next  to  each  other  at  pixel  spacings  that  scale  relative  to  the  along  track  spacing  of  the 
pulses  Ax  =  vpr  where  vp  is  the  platform  velocity  and  r  is  the  pulse  repetition  period.  The 
final  image  is  essentially  a  raster  scan  of  a  strip  of  the  sea  floor,  hence  the  name  ’’strip-map 
image” .  Synthetic  aperture  imaging  is  a  coherent  imaging  technique  that  exploits  the  extra 
information  available  in  the  phase  of  the  real  aperture  data.  We  adopt  the  Stop  and  Go 
model;  a  point  source  radiates  at  time  t  =  0,  a  spherical  wave  that  reaches  the  sampling 
points  after  different  time  intervals.  If  the  source  is  placed  at  (x0,  z0)  the  time  t(x,  x0,  z0)  at 
which  the  wave  arrives  at  the  sampling  point  (x,  z)  is: 


The  field  d  due  to  a  distribution  s(x,  z)  of  source  emitting  at  t  —  0  can  be  expressed  by 


where  d  is  the  Fourier  transform  (in  time)  of  the  signal  d.  SAS  measures 


SAS(x,  t)  =  d(x,  z  —  0,  t) 


along  the  sonar  path  (x,  z  =  0)  =  T. 

Thus,  SAS  imaging  is  formulated  as  a  linear  inverse  problem; 

Problem:  Reconstruct  s(x,z)  from  SAS  data  SAS (x,t). 

Among  a  number  of  algorithms  mm  and  reference  therein,  which  have  been  developed 


for  Problem  the  frequency  domain  u-k  method  based  on  Stolt’s  map  m  mm  is  the  most 
efficient  and  accurate  method.  As  will  be  discussed  in  Section  4  it  has  certain  limitations,  es¬ 


pecially  it  assumes  the  homogeneous  scattered  media.  The  proposed  method  can  incorporate 
inhomogeneous  media  and  has  additional  capabilities,  (see  Section  4). 

An  outline  of  the  paper  is  as  follows.  In  Section  2  we  describe  a  geometric  migration 
method  based  on  one-way  wave  equations  for  reconstructing  s.  A  noble  algorithm  using 
the  integration  of  the  one-way  wave  equations  based  on  the  operator-splitting  method  is 
developed  and  its  stability  and  complexity  are  analyzed  in  Section  3.  In  Section  3  we  list 
advantages  of  the  proposed  method  comparing  to  the  uj-k  method.  The  image  enhancement 
technique  based  on  the  BV-type  reguralization  is  discussed  in  Section  4.  In  Section  5  we 
present  a  test  against  real  field  data,  provided  by  the  Naval  Surface  Warfare  Center-Panama 
City,  Florida  and  a  comparison  with  the  u-k  method. 
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2  Geometric  Migration 


We  construct  an  approximating  solution  based  on  the  geometrical  migration  via  the  one-way 
wave  equations.  Let 

A(kx,  uj)  =  Tx±  SAS(a;,  t). 

Assume  the  plane  wave  extrapolation 


with 


D(kxi  kz,  u)  =  A(kx,  u>)exp(j(u>t  +  kxx  +  kzz)) 


w2  =  -r(^  +  a. 


Then  the  inverse  Fourier  transform  of  D 


d(x,z,t)  =  ^  3  J  D(kx,kz,ou)dkxdkzdcu 


4  d2d  d2d  d2d 
+ 


satisfies  the  wave  equation 

^  c2dt 2  dx 2  '  dz 2 

with  the  boundary  condition  at  z  =  0 

d(x,  0,  t)  =  SAS(a:,  t) 


and 

~  f)(] 

d{x,  z,T)  —  0  and  zi  T)  =  0. 

Wave  equation  based  migration  integrates  the  wave  equation  (3)  backward  in  time  to  obtain 
an  approximation  s  of  distribution  s  as; 

d(x,  z,  0)  =  s(x,  d). 

SAS  data  is  created  by  integrating  over  the  beam-width  of  the  sensor.  The  radiation 
pattern  of  any  dimension  (width  or  length)  of  an  aperture  has  an  angular  dependence  that  is 
referred  to  as  the  beam  pattern  of  the  aperture.  Beam  patterns  are  frequency  dependent  and 

Q 

have  beam-widths  given  by  the  3dB  response  of  their  main  lobes;  6  =  aw  where  D  is  the 
length  of  the  aperture  and  /  are  the  frequency  of  the  signal  that  the  aperture  is  transmitting 
or  receiving.  The  term  aw  is  a  constant  reflecting  the  main  lobe  widening  due  to  weighting 
of  the  aperture  illumination  function.  For  example  /  =  120kHz  and  D  =  0.04m  and  aw  =  1 
gives  9  =  17.9  degrees  Thus,  in  order  to  speed-up  the  wave  equation  based  algorithm  and 
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minimize  the  undesirable  effects  of  unwanted  reflections  we  use  the  (15  degree)  one-way  wave 
equation  based  on 


(4) 


kz  =  k\ 


t) 


where  k  =  2 u/c  and  we  assumed  \kx/k\  «  1.  In  time  domain  (4)  is  equivalently  written  as 

1  d2u 


(5) 

with 


4  d2u  2  d2u 

c2  d t2  c  dzdt  2  dx2 


u(t,x,  0)  =  SAS (x,t),  x  e  r. 

An  advantage  of  the  method  is  that  it  allows  one  to  have  a  specified  variable  wave  speed 
c  =  c(x,  z )  of  media.  The  corresponding  method  for  the  polar  and  cylindrical  geometry  is 
given  as 

Polar  coordinate 


Cylinder 


4  d2u  2  d2u  1 1  d2u 
c2  d t2  c  dud t  2  r  d62 

4  d2u  2  d2u  1  ,ld2u  d2u N 


c2  dt 2  c  dudt  2  r  dO2  dz2 

We  can  derive  the  wide  angle  one-way  wave  equation  based  on  the  rational  approximation 


(6) 

we  have 


i  i  kx\2  a(kx/k)2  , 

kz  =  k\  1-  [-P-]  ~k(  1 - 

A/  k  J  1  1  —  P{kx/k)2' 


kz{k  —  ~rk2)  =  k2  —  (a  +  /3)/c, 

K 


The  differential  form  is  given  by 

(7) 


4  <92w  <9  f2du  c 

c 2  dt2  +  dz  \  c  dt  ^2 


(92u  \  d2u 

dt)  =  (a  +  P)- 


dx2  J  v  dx 2 

With  a  —  .5,  (3  —  .25  and  a  =  .478,  (3  =  .376,  (7)  is  called  45  degree  and  65  degree 
approximation,  respectively. 


3  Migration  by  the  operator  splitting 

Q 

With  normalization  of  the  time  (t)  by  the  wave  speed  -  and  reverting  the  time,  (5)  is  written 


as 

f  Ut\ 

o 

o 

(S)  = 

d 

V  v‘ ) 

V  0  / 

u 


+ 


/  0  1  \ 
i  d2 

V  2  fa2  0  / 


u 


4 


So,  we  apply  the  time  splitting  on  [t,  t  +  At]  of  the  Lie- Trotter  form 


The  first  step  of  (9)  is  equivalent  to  the  shift  operation; 


{v(t  +  At,  x,  z )  =  v(t,  x,  z  —  At),  z  >  At 
d 

v(t  +  At,  x,  z)  =  — SAS  (x,  t  +  z),  z  G  [0,  At) 

The  second  step  of  (9)  is  the  one-D  wave  equation  in  x  and  is  well-posed.  In  fact,  let 

n  =  \—L,  L]  x  [0, 1]  and  H1,X(Q)  =  U  G  L2(fl)  :  e  L2(ty\.  Let  X1  =  x  L2Al) 

ox 

be  the  Hilbert  space  equipped  with 

C  r)n, 

\{u,v)\2Xl  =  /(|-|2  +  2|u|  )2)dxdz. 

Dehne  the  linear  operator  A\  on  X\  by 


Ai{u,i ;) 


1  d2u 

2  dx2 


dom(Ai)  =  {(u,u)  G  X1  :  v  G  Hl'x( ft),  —  G  L2(fi) 

ox1 


with  —  (±L,  z)  =  0} 
ox 

Then,  Ai  is  dissipative  and  skew-adjoint  on  Ad  and  thus  generates  a  strongly  continuous 
group  on  Xi.  Hence,  it  is  easy  to  show  that  if  (u,  v)  is  generated  by  (9)  then 

ft+At  g 

\(u,v)(t  +  At)\2Xl  <  |(«,u)(i)|x1  +  |— SASfA,  s)|2  dxds 


\(u,v)(T)\2Xi  <  /  |— SAS (x,t)\‘ 


Similarly,  we  can  argue  that  (8)  itself  is  well-posed,  i.e.,  if  we  dehne  the  operator  A  on  Ah 

by 

A{u,v)  =  (v,divXtZ(-  —  ,-v)) 

2  ox 


1  du  du 

dom(A)  =  {(u,v)  G  Ah  :  v  G  H1'x( ft,),  divXjZ(-  —  ,  —v)  G  L2{Vt)  with  v(x,z)  =  0,  —(±L,z) 

A  kJ ,1  U,L 
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then  A  is  dissipative  and  generates  a  contractive,  strongly  continuous  semigroup  on  Xx. 
We  fully  discretize  (9)  and  obtain 

Algorithm  I 

SAS+1  -  SAS!1 


+++,  =  vX  1  <  J  <  n  with  +0+1  = 


,n+ 1 


(10) 


At 


u 


n+l 


U 

-(/+ +  AW).  A1  =  - 


n+l 


U 


At 


-,  1  <  j  <  min(n,  M ) 


where  uX  and  vX  represents  the  value  of  u  and  v  at  the  grid-point  (iAx,jAz)  at  time  nAt, 


hi 


A  z 


respectively.  Here,  At  =  A z  and  c  =  — — ,  H  e  j^n+i.n+i  tri-diagonal  matrix  dehned 

Z— \<X/ 

by 

{Hu)i  =  ~(ui+i  -  2  Ui  +  2  <i  <  N, 


and  (Hu) i  —  ~(u2  —  Ui),  (iTu)jv+i  =  %+i  -  mat 


and  corresponds  to  the  central  difference  approximation  of 


d  2u 
dx2 


.  Also,  we  used  the  implicit 


Euler  scheme  to  integrate  the  second  step  (1-D  wave  equation  in  x).  That  is, 

,n+l  _  ~n 


u?+1  ~  U? 


1  r(  Hu),,  V" 


=  u 


n+l 


At  Ax2  At 

The  number  of  operations  at  the  n-th  time  step  of  (10)  is  of  order  0(N  min (n,M)).  M  is 
the  number  of  the  focusing  step  at  each  pixel  (i,j)  in  cross-range  direction  x  and  if  j  >  M , 
then  =  +,  Thus,  the  total  operation  is  of  order  0(MM2). 

For  the  wide  angle  equation  (7)  we  define 


F=--~f 3- 
cdt  2 


rt  d2u  ,  2d u 

at  and  v  = 


It  follows  from  (7)  that 


and 


2  OF 
c  dt 


Thus,  (7)  is  equivalent  to 


'0 


dx2 


,2  2  <9+  < d2u 

(3- 


■ c '  dt 2 
d 


dx2 


c  dt 


&F  (Pu 

dz  +  a  dx2 


2dF  dF 
c  dt 


d2u 


dz  a  dx2 


(11) 


2  du  . 

cM={v-F)  +  F 


2d,  n  d2u 

cdt^’~  )  = 
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With  v  —  v  —  F,  we  use  the  three  step  splitting: 


(12) 


(  dF_  dF_  _ 
~dt+~dl~ 

du 

~dt  = 
dv 

dt  = 


(  dF  d2F 


du 

dt 


=  v 


dv  d2v 
dt  =  '  OF 


dt  a  dx 2 


du  _  F 
m~F 


dv 

dt 


=  0 


If  (3  —  0  then  v  =  0,  F  =  v  and  thus  it  reduces  to  the  two-step  splitting  method  (9).  The 
first  equation  is  accompanied  by  the  boundary  condition 

F(t,  x,  0)  =  J^-SAS(i,  x)  —  v(t,  x ,  0). 

Each  step  of  (12)  is  a  well-posed  linear  system  as  shown  above  and  we  can  prove  that  (11) 
is  well-posed.  In  fact,  let  =  [— L,  L]  x  [0, 1]  and  define  the  linear  operator  on  A2  on 
X2  =  L2(Q)  x  H1,X(Q)  x  L2{ ft)  by 

.  , _ j  .  .  dF  d2u  n d2u 

A2(F,u,v)  =  (-—  +  a^-2,v  +  F,P  — ) 


with 


d 

dom(A2)  =  {-q-F  e  with  F(-,0)  =  0  and 

fp  r)i  i  c) 

—u  G  L\n)  with  -(±L,  z)  =  0,  -(u  +  F)  E  L2m- 


We  equip  A"2  with  norm 


\(F,u,v)\x2  =  JJ, \-^;u\2  +  ~\F\2  +  jj\v\z)  dxdz 


Then,  A-2  is  dissipative,  i.e., 

(A2(F,u,v),  ( F,u,v )) 


f  ,d2u  du  d  dF 

=  L(xfl{v  +  F)+diTF  +  F)-  &  F) dxdz 


| F(x,  1)|2  dx  <  0. 


l-L 


Since  range(A2 )  =  X2,  A2  generates  a  strongly  continuous,  contraction  semigroup  on  X2. 
Similarly,  we  have  the  energy  estimate 


(I  J^eoi2  +  ):iF(r)i2  +  ^\v(t)\*  dxdz  < 


| F(t,  x,  0)|2  dx  dt. 
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Algorithm  I  is  extended  to  integrate  (12)  as  follows; 

Algorithm  II 


QAQn+1  _  c  a  cjn 

iffi1!  =  Kj,  1  <3<n  with  FI ^  =  bAb-  M  bAb-  -  5?0 


u 


n+1 

■J 


=  (/  +  (5c2H)~\uF  +  AtF.^1),  F.71  =  -J  ' J  ,  1  <  j  <  min(n,  M) 


un+l  _  ^n+1 

i#1  =  (/  +  ac2H)~l(un+l  +  Aft;™),  <j+1  =  ',J  A ^  ,  1  <  j  <  min (n,  M). 


At 


That  is,  we  require  double  the  operations  for  the  integration  of  Algorithm  II. 


4  Advantages  of  the  proposed  methods 

The  frequency  domain  ui-k  method  based  on  Stolt’s  map  [7]  is  the  most  efficient  and  accurate 
method  for  the  homogeneous  media  due  to  the  efficiency  of  fast  Fourier  transform.  It  also 
assumes  a  rectilinear  sonar  path. 

We  can  use  our  proposed  algorithms  as  a  means  to  compensate  the  motion  of  sonar  path. 
That  is,  let  T  be  a  curved  sonar  path  and  To  is  a  reference  rectilinear  path  (z=0).  Then  we 
solve  (5)  or  (7)  on  the  domain  enclosed  by  the  boundaries  T  and  To  with  boundary  value 

u(t,x,z )  =  SAS (t,x),  (x,z)  G  T 

In  this  way  we  have  the  mapped-SAS  data  u(t,x,  0)  at  T0  and  then  apply  the  omega  —  k 
method  for  the  rectangular  domain  hi. 

Our  implementation  (10)  of  the  one-way  wave  equations  is  easily  adjusted  to  the  case  of 
layered  media  c  =  c(z)  by  varying  the  range  increments  A z. 

The  proposed  method  can  allow  to  localize  the  integration  on  sub-layered  regions  (assum¬ 
ing  the  homogeneous  media).  Also,  we  can  integrate  (5)  or  (7)  in  overlapped  sub-domains  in 
the  cross-range  (x)  direction  and  then  apply  the  superposition.  This  improves  the  efficiency 
of  the  proposed  algorithms. 

5  BV-type  Regularization  for  Enhancement  of  SAS  imag¬ 
ing 

SAS  imaging  s(x,  z )  may  be  altered  by  inhomogeneity  of  the  field,  sensor  noise  and  irreg¬ 
ularity  of  the  sonar  path  and  so  on.  We  use  the  image  enhancement  technique  based  on 
BV-type  reguralization  [5]. 


Enhancement  S  of  s  minimizes 


(12) 

where 

and 


r 

\S  —  s\2  dxdz -\- (3  /  <^(|  —  |2 + 

Jo  <JX 


3S 

dz 


dxdz 


(3  >  0  is  the  regularization  parameter 


Q{4>)  =  /  <^(|  VS1)2)  dxdz  defines  the  restoration  energy. 

Jn 

The  followings  summarize  our  findings  in  [5]  on  the  enhancement  based  on(12); 

•  tp(t2)  =  t2  corresponds  to  the  standard  Gaussian  filter  and  works  well  for  a  smooth 
image  s. 

•  <p(t2)  =  t  corresponds  to  the  BV  (nonlinear)  filter  and  restores  edges  and  flat  regions 
of  image  s  very  well.  But,  it  has  significant  stair-case  effects. 

•  In  order  to  deal  with  images  with  multi-scales  of  edges,  flat,  and  smooth  regions  we 
developed  an  algorithm  which  uses 
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s  e  1,00) 
1  s  e  [5, 1] 

—/=  s  £  (0;  <5) 

Vs 


v'(s)  =  < 

< 

It  is  based  on  the  scale  analysis  and  we  demonstrated  the  applicability  of  the  algorithm 


m 


•  The  necessary  and  sufficient  condition  of  (12)  is  given  by 

-f3V  ■  &'(\VS\2)VS)  +  S  =  s. 

An  efficient  algorithm  for  finding  S  based  on  the  fixed  point  iterate; 

-(3  V  •  {(p\\X7Sk\2)VSk+1)  +  Sfc+1  =  s 
is  developed  and  analyzed  in  |5)  and  is  used  in  our  test. 


6  A  Test 

The  algorithm  is  successfully  applied  to  real  data  that  are  available  to  us  via  the  Naval  Surface 
Warfare  Center  (NSWC)  and  shows  a  promising  capability.  A  full  capability  is  going  to  be 
tested  in  the  line  of  its  advantages  discussed  in  Section  4.  In  a  CRSC  tereport,  CRSC- 
TR09-12  at  http://www.ncsu.edu/crsc/reports/reports09.htm  we  show  the  raw  SAS  data, 
SAS  imaging  by  algorithm  (10),  and  the  image  enhanced  by  our  enhancement  algorithm. 


9 


References 


[1]  J.F.  Clarebout,  Coarse  grid  calculations  of  waves  in  inhomogeneous  media  with  appli¬ 
cations  to  delineation  of  complicated  seismic  structure,  Geophysics,  35  (1970),  407-418. 

[2]  M.N.  Guddati  and  A.H.  Heidari,  Migration  with  arbitrary  wide-angle  wave  equations, 
Geophysics,  70  (2005),  S61-S70. 

[3]  D.W.  Hawkins,  Synthetic  aperture  imaging  algorithms:  with  applications  to  wide  band¬ 
width  sonar,  Ph.D  thesis,  University  of  Canterbury,  1996. 

[4]  M.P.  Hayes  and  P.T.  Gough.  Broad-band  synthetic  aperture  sonar.  IEEE  Journal  of 
Oceanic  Engineering,  17  (1992),  80-94. 

[5]  K.  Ito  and  K.  Kunisch,  BV-type  Regularization  methods  for  convoluted  objects  with 
edge- flat-grey  scale,  Inverse  Probles  16  (2000),  909-928. 

[6]  M.  Soumekhi,  Fourier  Array  Imaging,  Prentice  Hall,  Englewood  Clifs,  NJ,  1994. 

[7]  R.H.  Stolt,  Migration  by  Fourier  transform,  Geophysics,  43  (1978),  23-48. 


10 


